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Abstract. A grand canonical Monte Carlo method for the simulation of a simple 
colloid-polymer mixture called the AO model will be described. The phase separa- 
tion known to occur in this model is driven by entropy. The phase diagram of the 
unmixing transition, the surface tension and the critical point will be determined. 



1.1 Introduction 

Mixtures of particles are all around us. If you pour oil and water together you 
create a mixture. Another example is a solution containing colloids and poly- 
mers. Mixtures will sometimes phase-separate or unmix. When this happens 
particles of the same kind cluster together. A famous example is a mixture of 
oil and water: after unmixing, the lower part of the container contains mostly 
water and the upper part mostly oil with an interface in between. 

There are different reasons for a mixture to unmix. For instance, the 
unmixed state might carry a lower energy. As an example consider a mixture 
of A and B particles interacting with the following pair potentials: 

UAA(r) = u BB {r) = 0, u AB (r) = -, (1.1) 

r 

where r is the distance between two particles and e some positive constant. 
In this mixture, A particles do not feel each other, and neither do B particles, 
but when A and B particles are close together there is an energy penalty. At 
low temperature, the system tries to minimize the energy and can only do so 
by moving the A and B particles as far apart as possible. The system thus 
unmixes. 

There is another mechanism that can induce unmixing. This mechanism 
has its origin in entropy and not in energy The unmixing of colloid-polymer 
mixtures for example is driven by entropy. To show that this must be the 
case, it is instructive to briefly consider the nature of the interactions in a 
typical colloid-polymer mixture. To a crude approximation the colloids be- 
have as hard spheres. The polymer interactions are more complicated. A real 
polymer consists of a chain of bonded monomers. In principle all monomers 
should be taken into account explicitly. However, under certain conditions 
it is reasonable to describe the polymer chain with only the coordinate of 
its center of mass and some effective radius called the radius of gyration. In 
1954 Asakura and Oosawa [1] proposed a simple model for colloid-polymer 
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Fig. 1.1. Colloid and polymer at minimum separation (here and throughout this 
text C=colloid and P=polymer). The distance between the centers of mass equals 
R c +Rp. This is the minimum allowed separation: any smaller separation is punished 
with infinite energy. The colloid thus carries a spherical depletion zone with volume 
Vs around its center of mass which cannot contain any polymer centers. 



mixtures based on precisely these approximations. In this model colloids and 
polymers are treated as spheres with respective radii R c and R p . Hard sphere 
interactions are assumed between colloid-colloid (cc) and colloid-polymer (cp) 
pairs while polymer-polymer (pp) pairs can interpenetrate freely. This leads 
to the following pair potentials: 



u pp (r) = 0, 

with r again the distance between two particles. The above equations define 
what is nowadays called the AO model. Note that in 1976 the same model 
was proposed again and independently by Vrij [2]. 

According to Ea. (|1.2|) the energy of a valid AO mixture is always zero. 
Therefore, if unmixing occurs (and it does) it cannot be explained in terms 
of the energy argument used before. The reason for unmixing is a little more 
subtle, see Fig ll.ll and Fig ll.2l Shown in Fig ll.ll is one colloidal particle just 
touching one polymer. When the particles touch, the distance between the 
centers of mass equals R c + R p . This is the minimum allowed separation 
between the particles because any smaller separation carries an infinite en- 
ergy penalty. In the AO model every colloid is thus surrounded by a region 
which cannot contain any polymer centers of mass. This region is called the 
depletion zone. In three dimensions the volume of the depletion zone equals: 



Consider now Fig ll.2l which shows a container of volume V holding zero, 
one and two colloidal particles. Shown at the top is a polymer which we 
want to insert into the container. If the container is empty we can place the 
polymer anywhere so the free volume / available to the polymer is simply 




(1.2) 



Vs = —{Rc + R P f. 



(1.3) 
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Fig. 1.2. Free volume / available to a single polymer in a container of volume V 
holding zero, one and two colloidal particles. See text for details. 



/ = V CFig |1.2l !V). If the container already holds one colloid the free volume 
decreases to f — V — Vg fFig ll.2B ). The physics becomes interesting when 
the container holds two or more colloids. Two cases can now be distinguished. 
In the first case both colloids are well separated and the free volume equals 
/ = V — 2Vs CFig |1.2Q . In the second case the colloids are so close together 
that their depletion zones overlap and the free volume increases / > V — 
2Vs (Fig ll.2D ). This immediately has physical consequences: by clustering 
together the colloids can increase the free volume available to the polymers 
and hence the entropy of the polymers. Under certain conditions the gain in 
entropy is sufficient to drive unmixing. 

The AO model thus provides a convenient framework in which to study 
entropy driven unmixing. It has sparked much theoretical work and many 
simulations [3-7]. Despite its apparent simplicity, predicting the phase be- 
havior of the AO model is no easy task. In this paper an efficient grand 
canonical Monte Carlo (MC) method that can be used to simulate the AO 
model will be described. The use of the grand canonical ensemble allows one 
to bypass certain problems encountered for example in the canonical and 
Gibbs ensembles. In particular, one can accurately obtain the surface tension 
and also perform simulations close to the critical point. The method will be 
used to calculate the phase diagram and the surface tension of the interface. 
We will also present an accurate determination of the critical point using 
finite-size scaling. 



1.2 Grand Canonical Monte Carlo 

Imagine an AO mixture contained in some volume V at inverse tempera- 
ture j3 = l/(fcsT). In the grand canonical ensemble V and (3 are fixed but 
the number of particles inside V is allowed to fluctuate. The probability of 
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Fig. 1.3. Problem encountered in a standard grand canonical MC simulation of the 
AO model. When the number of polymers in the volume is substantial it becomes 
difficult to insert additional colloids. Inserting more polymers remains easy though. 



observing a volume containing N c colloids and N p polymers is given by the 
grand canonical distribution: 

P = Cz?°z?>e- flE , (1.4) 

where C is a normalization constant, E the energy given by Eq. l|1.2|l and 
{z c , z p } the fugacities of the colloids and polymers, respectively. 

In a grand canonical MC simulation of the AO model one would like to 
generate configurations of colloids and polymers that sample Eq. i|1.4fl . In the 
standard approach this is done by attempting to insert a single particle into V 
at a random location, or remove a single (randomly selected) particle from V. 
Insertion and removal are usually attempted with equal probability and ac- 
cepted with probabilities that depend on the energy change of the attempted 
move and on the fugacities. The standard approach is however not efficient 
for the AO model as Fig ll.3l illustrates. The figure shows a volume containing 
a substantial number of polymers and some colloids. According to Ea. l|1.2fl 
polymers do not interact with each other so they are happy to overlap. This 
is what generally will happen as the figure shows. Unfortunately, it is nearly 
impossible to insert an additional colloid into this system: no matter where 
the colloid is placed it will likely overlap with at least some of the polymers, 
and colloid-polymer overlaps are forbidden by Eq. l|1.2|) . A standard grand 
canonical MC simulation of the AO model is thus characterized by a very 
low acceptance rate of colloid insertions. 

The problem illustrated in Fig ll.3l is typical of asymmetric binary mix- 
tures of which the AO model is an example. The standard grand canonical 
MC algorithm does not deal well with such mixtures, essentially because it 
moves only one particle at a time. A MC move capable of removing entire 
clusters of polymers would be much more efficient. By using such a cluster 
move the formation of "holes" in the "sea" of polymers is enhanced. If the 
holes are large enough to contain a colloid, the acceptance rate of colloid 
insertions will increase. The MC move used in this work to simulate the AO 
model is aimed at doing precisely that. 
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Fig. 1.4. Inserting an additional colloid into an AO mixture. (A) The starting 
configuration. (B) The colloid is inserted at a random location inside the volume. 
(C) All polymers that overlap with the colloid are removed. 



1.3 Cluster Moves 



Consider now Fig ll.4t \ which again shows an AO mixture of colloids and 
polymers. This time we are less shy and simply insert an additional colloid 
at some randomly selected location inside V. The resulting configuration is 
shown in Fig ll.4b . with the colloid inserted in the upper right corner of the 
box. Note that the configuration in Fig ll.4b is not a valid AO configuration 
because the colloid overlaps with four polymers. To fix this the overlapping 
polymers are simply removed, leading to the configuration shown in Fig ll.4C . 

The MC move of Fig ll.4l is of course not sufficient to carry out a simu- 
lation. We also need a reverse MC move capable of bringing us back from 
the configuration of Fig ll.4C to the starting configuration of Fig ll.4t \.. In 
constructing the reverse move, the key idea is that for every colloid you take 
out, a number of polymers must be inserted into the empty depletion zone 
left behind by the colloid. But exactly how many polymers? 

To answer this question it is instructive to consider an AO model con- 
taining only polymers. In this case the AO model reduces to a very simple 
system, namely that of an ideal gas (remember that the polymers do not in- 
teract with each other). According to basic statistical mechanics, the average 
density of an ideal gas is equal to its fugacity. Likewise, for an AO model free 
of colloids, the average polymer density equals z p . If we attempt to insert 
a colloid into the pure polymer system it will on average overlap with z p Vs 
polymers with fluctuations in the average that are Poisson like, i.e. of order 
yj ZpVg. Because of these fluctuations, the number of polymers n p that we 
must insert for every colloid that we take out cannot simply be a constant. 
Instead, n p must be a random variable drawn from some probability distri- 
bution. One choice that works well (see Section fOjl is to draw it uniformly 
from the interval n p G [0, to) (so including but excluding to) and to an 
integer given by: 

to = 1 + max 1, int (z p Vs + ot\J ZpVs^j , (1-5) 
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Fig. 1.5. MC move used to remove a colloid from an AO mixture. (A) The colloid 
to be removed is chosen randomly from among those present. (B) The colloid is 
removed and n p random locations are selected inside the empty depletion zone left 
behind. The number n v is uniformly drawn from the interval [0, m). In this example 
n v — 4. (C) Polymers are placed onto the random locations. 



with a a positive constant of order unity you are free to choose (a rs 2.0 
usually gives good results). Recall that V$ is the volume of the depletion zone 
given by Eq. (|1.3|l . 

The reverse move can now be constructed and is illustrated in Fig ll .51 
First, we randomly select a colloid (Fig ll. 51 ^. The colloid is removed and n p 
random locations inside the depletion zone of this colloid are selected, with 
n p uniformly drawn from the interval [0, m) (Fig ll. 5B ). Finally, polymers are 
placed onto these random locations resulting in the configuration shown in 
Fig[T3p. 

1.4 Detailed Balance 

We will now derive the acceptance rates of the MC moves illustrated in 
Fig ll. 41 and Fig ll. 51 The acceptance rates must be constructed such that 
detailed balance is obeyed. This is vital because the algorithm must sample 
the grand canonical distribution of Ea. l|1.4[) . Recall that the condition of 
detailed balance demands that [8-10]: 

P(a)g(a r)A(a -f r) = P{t) 9 {t -> <j)A(t -> a), (1.6) 

with P(cr) the Boltzmann probability of state a, g(a — > r) the probability 
that the MC scheme generates state r from state a, and A(a — > r) the 
probability of accepting the new state t. 

The derivation that follows is based on two assumptions: 

1. The integer m of Eq. l|1.5|) is assumed a constant parameter of the algo- 
rithm. It must be set once at the start of the simulation and it may not 
be changed during the course of the simulation. 

2. It is assumed that the insertion of a colloid and the removal of a colloid 
are attempted with equal probability. 



1 Entropy Driven Phase Separation 7 



In an implementation it is important that the above conditions are met. If 
they are not the acceptance rates to be derived next may yield wrong results! 

1.4.1 Colloid Removal 

The acceptance rate for the removal of a colloidal particle is derived first. 
This is the MC move shown in Fig ll.5l Assume that we start in a state a 
containing N c colloids, N p polymers and energy E a . It is convenient to label 
the state as cr(N c , N p , E a ). The energy E a will of course be zero but for the 
derivation it is convenient to write it down explicitly. After removing the 
colloid, n p polymers are inserted into the depletion zone left behind. The 
state r that we end up in can thus be labeled as t{N c — 1, N p + n p , E T ). Note 
that E T need not be zero: if a second colloid happens to be very close to the 
colloid that was removed, it could now overlap with one or more of the n p 
polymers that were just inserted. 

To enforce detailed balance the probabilities that appear in Ea. (|1.6|l must 
be written down. We begin with the easy ones P(a) and -P(r). These are 
simply given by the grand canonical distribution of Eq. I|1.4I) : 



The probabilities g(a — > r) and g(r — ► a) are more complicated. If we 
are in state a the probability of ending up in state r is given by the product 
of the probabilities of the individual steps that were taken in the MC move. 
After close inspection of Fig ll.5l the following steps can be identified: 

1. Selecting a colloid at random. Since state a contains N c colloids the 
probability of choosing one particular colloid is l/N c . 

2. Selecting n p . Since n p is drawn uniformly from the interval [0, m) the 
probability of this step is l/m. 

3. Selecting n p random locations inside a volume Vs. The probability of 
choosing one particular location equals 1/Vs (see Appendix 11.9(1 . Since 
n p locations are selected we must raise this probability to the power 
n p . Additionally, we pick up a factorial counting the number of ways in 
which the locations can be selected. The total probability of this step 
thus becomes {n p )\/Vg p . 

The probability g(a — * r) is thus the product of the above three terms. 

Finally, we derive g(r — > a). This is the probability that the reverse MC 
move brings us back from state r to state a. This move is illustrated in Fig ll.4l 
and upon inspecting it we observe that state a is regenerated if a colloid is 
placed inside V at precisely the same location it was just removed from. This 
carries a probability 1/V and so we find g(r — > a) = 1/V. 

Substitution of the above terms into Ea. (|1.6|> and using the Metropo- 
lis choice [8] we find that the removal of a colloid must be accepted with 
probability: 



P(r) = Czl 



N c -1 N p +n p -0E r 



■p 



(1.7) 



A(N C — ► N c — 1) = min 1, 



mN c {ZpVsY^ 
z c V (n p )\ 



(1.8) 



8 Richard Vink 



Note that moves leading to configurations where E T is not zero are automat- 
ically rejected by the above equation. 



1.4.2 Colloid Insertion 

To make the algorithm complete the acceptance rate for the insertion of a 
colloid must still be derived, i.e. the MC move of Fisr ll.41 We can follow 
the same reasoning as before with perhaps one subtlety. We start again in 
a state a containing N c colloids, N p polymers and energy E a , so with label 
a(N c , N p , E a ). All polymers (say n p of them) that overlap with the inserted 
colloid are removed, so the state r that we end up in can be labeled as 
t(N c + 1,N P — n p ,E T ). Again, the energy E T need not be zero: if we are 
unlucky the inserted colloid overlaps with one or more of the N c colloids 
already present. According to Eq. l|1.2|) such overlaps also carry infinite energy 
To construct detailed balance we first write down P(a) and P(r): 

P{a) = Cz^z^e-P E % P{t) = Cz^z^-^e^ . (1.9) 

Next, we derive g(a — > r) which is the probability that the MC move of 
Fig 1 1.41 generates state r starting in state a. This involves only the random 
selection of a location inside V so we obtain: g(a — ► t) = 1/V. 

Finally, we consider g(r — > a). This is the probability that the reverse 
MC move of Fier ll.51 transforms state r back into state a. This will happen 
only if the following conditions are met: 

1 . The newly inserted colloid is removed again. Since state r contains N c + 1 
colloids the probability of this step equals 1/{N C + 1). 

2. Precisely the same number n p of polymers that were removed are inserted 
again. This step requires some care. Assume that we chose m rather small 
in Eq. (|1.5() . In that case the forward move of Fig ll.4l might remove more 
polymers than the reverse move of Fig 1 1.51 can possibly put back. This 
will happen if n p > m. Therefore, the probability p that exactly the same 
number of polymers are inserted must be written as: 

_ f if n p > m , . 

[1 /m otherwise. 

3. Selecting the same n p coordinates inside V$. As was explained before this 
probability equals (n p )\/V^ p . 

The probability g(r — ► a) is thus the product of the above three terms. 

Substitution into Eq. I|1.6fl and using the Metropolis choice [8] we find that 
the insertion of a colloid must be accepted with probability: 

A{N C - JV C + 1) = (1.11) 

| if n p > m 



mm 



otherwise. 



m(AT c + l) (z p V 5 ) n P 

Again, moves leading to configurations where E T is not zero are automatically 
rejected by the above equation. 



1 Entropy Driven Phase Separation 9 




Fig. 1.6. Even low values of m may yield high polymer densities. (A) The starting 
configuration showing a colloid with its depletion zone. The colloid is removed 
and one polymer is inserted at a random location inside the depletion zone (cross) 
resulting in configuration (B). (C) Another colloid is inserted close to the polymer. 
(D) The colloid is removed again with one polymer inserted in its depletion zone 
at the cross resulting in two overlapping polymers. 



1.5 Ergodicity 

Now that the MC scheme obeys detailed balance we need to check if it is 
ergodic. An algorithm is ergodic if every point in phase space can be reached 
with finite probability. It is easy to see that the algorithm is ergodic as far 
as the colloids are concerned. For every inserted colloid, a random location 
inside V was selected so the entire volume is correctly sampled. Similarly, 
polymers also sample the entire volume but do so indirectly, namely via the 
removal of colloids. Whenever a colloid is removed, a number of polymers are 
inserted in the depletion zone. Since the colloids sample the entire volume, 
so then do the polymers. 

One might argue that the choice of to in Ea. (|1.5fl violates ergodicity. If 
at most to — 1 polymers are inserted for every colloid, the polymer density 
will always be less than (m — 1)/Vs as a result. While this objection sounds 
reasonable it is in fact unjustified as Fig 1 1.61 shows. 

Fig |1.6l shows an example simulation of the AO model using the cluster 
algorithm just described with a low value of m, namely to = 2. This means 
that for every removed colloid, at most one polymer is inserted. The starting 
configuration is a volume containing one colloid, see Fig ll .6t A. The edge of the 
depletion zone of the colloid is also drawn. In Fig ll.6b the colloid is removed 
and one polymer is placed very close to the edge of the depletion zone (but 
still inside of it, of course). The simulation is continued in Fig ll.fiC where 
another colloid is placed into the system. This colloid is placed close to the 
polymer already present, but not too close so the polymer survives. Finally, 
in Fig ll.BL ) the colloid is removed again and one polymer is placed inside the 
depletion zone of this colloid, close to the other polymer. The configuration 
in Fig ll.GD now shows two polymers practically stacked on top of each other, 
even though the algorithm was run with to = 2. This removes the above 
objection. 

Running the algorithm with such a low value of m is not recommended 
though. The whole point of this discussion is to show that m does not influ- 
ence the correctness of the algorithm, only its efficiency. 



10 Richard Vink 



1.6 Early Rejection Scheme 

In Section fl.3l it was argued that the number of polymers n p that one needs to 
insert for every colloid that is removed cannot simply be a constant. In fact, it 
was shown that n p is a random variable described by a Poisson distribution. 
Yet, in the subsequent description of the algorithm it was decided to draw 
n p uniformly, and not from a Poisson distribution. One might wonder if the 
efficiency of the algorithm is not seriously impeded by this choice. The answer 
is it is not, provided one implements the so-called early rejection scheme. 

In many MC simulations, the usual approach is to make a change to 
the system, calculate the involved energy change, select a random number 
between zero and one, compare this random number to the acceptance rate 
and finally accept or reject the move. In some cases, particularly in systems 
where the interactions are hard sphere like, one can do much better. 

Consider for example the removal of a colloidal particle displayed in 
Fig ll.5l and the associated acceptance rate Eq. l|1.8|) . Most of the quantities 
appearing in Ea. (|1.8() are already known at the start of the move. For exam- 
ple N c , V$, E a , to, [3 and the fugacities. In fact, the only unknowns are n p 
and E T . We also know that if the move is ever going to be accepted E T will 
be zero at the end. So the only remaining unknown is n p which, if you recall, 
is a random number drawn from the interval n p 6 [0, m) . The early rejection 
scheme proceeds as follows: 

1. Select a random number r between zero and one. 

2. Select n p uniformly from the interval n p £ [0, m). 

3. Calculate the acceptance rate A(N C — > N c — 1) using Ea. l|1.8|l with the 
above value for n p and assuming that E T is zero. 

4. If r > A(N C —> N c — 1) reject the move immediately otherwise proceed 
to the next step. 

5. Remove the colloid and insert the n p polymers. If any of the inserted 
polymers produce overlap with other colloids reject the move, otherwise 
the move is accepted. 

The most CPU time consuming step in the above scheme is step five. How- 
ever, this step is only performed for those values of n p that are reasonable. 
Unreasonable values for n p were already filtered out in step four at the cost of 
only a few multiplications and selecting two random numbers. In this scheme 
it therefore does not matter so much what distribution n p is drawn from. 
Needless to say, the early rejection scheme is highly recommended. To speed 
up the determination of overlap the link-cell method should also be used. 

1.7 Application 

In this section the grand canonical cluster scheme will be used to study bulk 
phase separation in the AO model. Since the AO energy is either zero or 
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Fig. 1.7. Probability P(r] c ) of observing a colloid packing fraction r\ c for an AO 
mixture with {q = 0.8; z c = 87.2; rf v — 1.0} in a simulation of dimensions L x = 
L y — 16.7 and L z = 33.4. The probability is not normalized. The inset shows the 
logarithm of the probability distribution. 



infinity the temperature plays no role so we simply put /3 = 1 in the accep- 
tance rates. The phase behavior of the AO model is thus fixed by the colloid 
to polymer size ratio q = R p / R c and the fugacities {z c , z p }. We consider here 
a size ratio q = 0.8 and put R c = 1 to set the length scale. The simulations 
are performed in a box with edges L x x L y x L z and using periodic boundary 
conditions. Following convention we define 77^ = z p (47r/3)i?p known as the 
polymer reservoir packing fraction. Note that 77^ has no relation to the num- 
ber of polymers actually in the system. It is just a different way of expressing 
the polymer fugacity z p . 

In a naive implementation of the scheme one sets the fugacities {z c ,r]p} 
and starts the simulation. As the simulation proceeds colloids and polymers 
will enter and leave the box. The crucial quantity to measure is P(j] c ) de- 
fined as the probability of observing a box with colloid packing fraction 
rj c = (Att/3)R^.N c /V . During the simulation one thus maintains a histogram 
counting how often a certain colloid packing fraction has occurred. 

If phase separation occurs P{i] c ) is bimodel. An example distribution is 
shown in Fig ll.7l The peak at low r\ c corresponds to the colloid vapor phase 
(V), the peak at high rj c to the colloid liquid phase (L) and the region in 
between is the phase-separated regime. The distribution in Fig ll.7l is at coex- 
istence which means that the area under both peaks is equal. At coexistence, 
the simulation spends equal time in both phases on average. Also shown is 
the logarithm of P(ry c ). The physical significance of this curve is its relation 
to the free energy. The height of the barrier marked AF corresponds to the 
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Fig. 1.8. Average polymer packing fraction r\ v as a function of the colloid packing 
fraction 77c for an AO mixture at coexistence with q = 0.8 and rfp — 1.0. 



free energy barrier separating the phases. As was shown by Binder [11] this 
barrier is related to the surface tension via 7 = AF/(2A) with A the area of 
the interface. 

It is also interesting to measure the average polymer packing fraction rj p = 
(4ir/3)RpNp/V as a function of r) c . This result is shown in Fig ll.8l In the pure 
polymer phase (r/ c = 0) we expect rj p = rfp because the polymers then mimic 
the ideal gas. This is precisely what Fig ll.81 shows. As the colloid packing 
fraction increases the polymer packing fraction in the system decreases to 
zero. 

In the AO model rf is the control parameter, much like temperature is 
for fluid-vapor transitions. To obtain the phase diagram one simply has to 
measure P(j] c ) at coexistence for a number of different rfp. The problem in a 
simulation is finding the value of z c that yields coexistence for the chosen rfp 
of interest. Additionally, if the barrier AF is high, it will be difficult to sample 
P{r]c) in the region between the peaks. Fortunately, an array of techniques 
is available to overcome these problems [12]. The results in this work were 
obtained using a new technique called Successive Umbrella Sampling [13]. 

For each P(rj c ) at coexistence one reads off the colloid packing fraction 
of the vapor phase r)Y and of the liquid phase r/^ and plots the two points 
(VciVp) an d (VciVp) m a graph. This yields the phase diagram in reservoir 
representation shown in the inset of Fig ll.9l By using Fig ll.8l we can convert 
the reservoir representation into the experimentally more relevant {rj c , rj p } or 
system representation also shown in Fig ll.9l For every P(r) c ) one also obtains 
a value for the surface tension using the method of Binder. The results of 
this procedure are shown in Fig ll.101 in two different representations. 
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Fig. 1.9. Phase diagram of the AO model with q = 0.8 in system representation. 
Crosses were obtained using box dimensions L x — L y = 16.7 and L z = 33.4; 
open circles were obtained in a smaller box with dimensions L x 
L z — 26.5. The inset shows the phase diagram in reservoir representation. 



Finally, we determine the critical polymer fugacity denned as the value 
of rfp above which phase separation begins to take place. From the inset 
of Fig ll.9l we see that the critical fugacity is around rf v ss 0.75. We have 
performed a finite size scaling analysis [14] by measuring the cumulant ratio: 

as a function of rf v close to the critical point for different system sizes. The 
results are shown in Fig ll.lll The critical fugacity is at the intersection of 
the lines from which we obtain rf v cr — 0.766 ± 0.002. 



1.8 Conclusions 



The grand canonical cluster scheme described in this chapter is very successful 
at modeling phase separation in AO mixtures. In fact, at the time of writing, 
it is unsurpassed in speed and accuracy by other simulation methods [15]. 
However, the method does have its limitations. If the packing fraction of the 
colloids is high (say 0.40 and above) the algorithm is no longer efficient. In 
that case the insertion of colloidal particles fails, not because of overlap with 
polymers, but because of overlap with other colloids. This problem is well 
known in hard sphere simulations. 

The algorithm could be improved for systems where the polymer-polymer 
interaction is not zero. In these cases the random insertion of n p polymers 
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Fig. 1.10. Reduced surface tension 7* = 4J?^7 f° r an AO mixture with q = 0.8 as 
a function of the difference in the colloid packing fractions in the coexisting liquid 
and vapor phases. The box dimensions were L x = L y — 16.7 and L z — 33.4. The 
inset shows 7* as a function of rf v . 
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Fig. 1.11. Cumulant ratio M given by Ea. lfl.12p as a function of for an AO 
mixture with q = 0.8 for various system sizes. The simulations were performed in 
a cubic box with edge L as indicated. From the intercept (vertical line) we obtain 
for the critical polymer fugacity r[ v ^ cr = 0.766 ± 0.002. 
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into the depletion zone of a colloid may no longer be efficient. Fortunately, 
the algorithm is easily adapted to include smarter insertion moves such as 
configurational bias [16], recoil growth [17], wormhole MC [18] and perhaps 
PERM [19]. Work along these lines is in progress. 

Note also that the presented method is general. The derivation of detailed 
balance for example can easily be modified to an AO model confined between 
walls or to systems interacting with smooth potentials. This would be the 
subject of further work. 
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1.9 Appendix: Random Points 

Here we argue that the selection of a random point inside a volume V carries 
a probability 1/V. The question is most easily answered by assuming that 
the volume is made up of many tiny cells, each with a volume r. The total 
number of cells in the volume is thus equal to n = V/T and the probability 
of picking one of them is l/n = r/V. 

The difficulty lies in choosing F. In statistical mechanics it is assumed that 
the smallest element of phase space has a volume determined by the thermal 
wavelength A leading to F = A 3 . This explains why the acceptance rates 
of grand canonical MC schemes often contain the thermal wavelength [9, 
10]. Unfortunately, in a MC simulation this is not very useful. Since the 
thermal wavelength depends on temperature, particle mass and even the 
Planck constant, this choice suggests that a simulation of for example hard 
spheres is temperature and mass dependent. 

The key observation is that the physics of the system is indifferent to the 
choice of F. To see this consider Ea. i|1.8fl which is the acceptance rate for 
the removal of a colloidal particle. Strictly speaking, the volumes V and Vg 
appearing in Ea. (|1.8(l must be replaced by V/T and Vg/T, respectively. These 
additional factors of F are, however, readily absorbed into the fugacities z c 
and z p . The fugacity is given by z — exp(/3/x) with fj, the chemical potential, 
so a rescaling of the fugacity merely shifts by a constant. This constant 
has no physical consequence because it in turn shifts the Hamiltonian by a 
constant. Thus, F has no physical consequence and may therefore be set to 
unity which was done throughout this text. 
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